1 N-PHASE INTERFACE TRACKING METHOD UTILIZING 

2 UNIQUE ENUMERATION OF MICROGRID CELLS 

3 

4 TECHNICAL FIELD 

5 

6 This invention relates generally to methods for the simulation of multiphase 

7 fluid flows, and more particularly, to those methods which track interfaces 

8 between immiscible fluids. 
9 

10 BACKGROUND OF THE INVENTION 

11 

12 The study of interfacial flow is a broad topic of interest in many different 

13 research disciplines. Physicists, biologists, engineers and other scientists all 

14 share a stake in accurate representation of interfacial position. Whether 

1 5 pipeline flow of crude oil mixtures, cellular disruption or the development of 

16 galaxies is to be modeled, interfacial movement intimately influences 

1 7 calculations by accounting for a redistribution of local density functions. In 

18 general, the problems of interest have as a common element the requirement 

19 to resolve through numerical simulation the complex flow patterns that result 

20 from the flow of immiscible (or semi-immiscible) fluids. 
21 

22 Published works dating from the 1960*s until now relate an ever-improving 

23 understanding of numerical algorithms that allow for an accurate description 

24 of interface position. These algorithms are generally categorized as either 

25 surface methods or volume methods. However, there Is still no generalized 

26 method presently developed that allows for the numerical simulation of the 

27 dynamic movement of three or more fluid materials and their interfaces. 
28 

29 A reference field or function that moves with an interface typically 

30 characterizes surface methods. Sometimes the reference is a mass-less 

31 fixed particle along the interface (Daly, 1969; Takizawa et. al, 1992). Other 

32 times a level set function which tracks the shortest distance to an interface 
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1 from a fixed point is offered (Sussman et. al, 1994; Osher and Sethian, 1988; 

2 Sethian, 1996). Further references indicate height functions are implemented 

3 where a reference line or plane is chosen in the computational domain 

4 (Nichols and Hirt, 1973). All of these methods give accurate descriptions of 

5 two fluid, single interface problems that do not involve folding, breaking or 

6 merging interfaces. The biggest advantage of surface methods is that an 

7 interface position is explicitly known for all time. 
8 

9 Volume methods, however, build a reference within the fluids under 

10 evaluation. Typically a method will discretely identify different materials on 

1 1 either side of an interface. Mixed cells then indicate the general location of an 

12 interface. An exact location is never discretely known, and volume methods 

13 are characterized by a reconstruction step whereby an approximate interface 

14 is built from local data consisting of volume and area fractions. 
15 

16 The volume methods may be further divided into 2 sub-cases to include 

1 7 particle methods and scalar advection methods. Particle methods include the 

18 marker and cell (MAC) method (Harlow and Welch, 1965; Daly, 1967) and the 

19 particle in cell (PIC) method (Harlow et.al, 1976). These algorithms employ 

20 the idea of mass-less particles to identify a particular material and then track 

21 their movements through a static grid based on local velocity conditions. The 

22 volume of fluid (VOF) method (Hirt and Nichols, 1981 ; Ashgriz and Poo, 1991; 

23 Lafaurie et.al. 1994) and various line techniques including SLIC (Noh and 

24 Woodward, 1976), PLIC (Youngs, 1982) and FLAIR (Ashgriz and Poo, 1991) 

25 implement several unidirectional sweeps to predict the change in volume 

26 fractions of cells at each time step based on a local resolution of the scalar 

27 transport equation. There are many combination algorithms presently used 

28 which implement both VOF and line techniques to capture even greater 

29 interface detail (Rider and Kothe, 1998; Gueyffier et al.. 1999; Scardovelli 

30 et al., 2000. 2002). Volume methods offer good interfacial descriptions of 

31 complicated fluid geometries in both two and three dimensions and allow for 
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interface folding, breaking, and merging. IHowever these methods are 
restricted to cases involving only two materials. 

Traditional VOF methods define a concentration function, C, to denote 
materials. Typically, 



where values between 0 and 1 represent mixed cells. 

C is then transported by the velocity field u via the scalar transport equation. 



And finally, the average (or cell centered) values of density, p , and 
viscosity, ^ , are interpolated as: p = C/?, + (1 - C)p^ and // = C//, + (1 - Qju^ 
where the subscripts 1 and 2 refer to materials 1 and 2. In this interpretation, 
C, although called a concentration function, acts as the fluid fraction of 
material 1 or 2 In any given cell. 

Methods have been employed previously which have attempted to track two 
interfaces or three materials in a computational fluid simulation. These 
methods (Eulerian-Eulerian) operate by using control volumes which may 
contain at most three materials and two unique interfaces. As discussed 
above, at least a single set of conservation equations is still required to 
resolve the flow field dynamics, i.e., mass conservation equation and two or 
three momentum conservation equations dependent on the dimensionality of 
the problem (2D or 3D), or a complete set of consen/ation equations is solved 
for each material. However, no one to date has addressed the issue of 
uniquely Identifying and tracking, with a single or multiple concentration 
function(s) Ci, three or more materials. The issue always distills down to the 
unique tracking and resolving of multiple interfaces in a computational cell. In 
addition, when multiple equation sets are used as possibly in Eulerian- 
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1 Eulerian methods, the numerical problem becomes intractable when folding, 

2 breaking, and reforming interfaces exist due to the need to precisely define 

3 the interface shape Inorder to apply the constitutive relationships between 

4 materials in this formulation. The above described computational fluid 

5 dynamic methods fail to provide a method which can readily handle explicit or 

6 implicit tracking of the interface between three or more generally immiscible 

7 fluid materials. The present invention provides an efficient and tractable 

8 method which overcomes the deficiency of current methods in tracking 

9 multiple fluid materials and their interfaces. 
10 

11 SUMMARY OF THE INVENTION 

12 

13 A method for tracking N materials and their interfaces in a computational 

14 domain is disclosed. A macrogrid including control volumes Is created on a 

15 computational domain in which N materials and their interfaces are to be 

16 tracked. A microgrid including microgrid cells Is overlaid upon the macrogrid 

17 with each of the microgrid cells being coupled to a control volume. The 

18 macrogrid and control volumes are initialized with initial and boundary 

19 conditions representative of the problem to be solved. A unique identifier 

20 number is assigned to each of the N materials and to the microgrid cells 

21 containing those N materials. Volume fractions for the N-materials in the 

22 control volumes can then be calculated. Similarly, average density and 

23 average viscosity for each of the control volumes may be discretely 

24 determined. 
25 

26 Equations of fluid motion (the Navier Stokes equations) may then be solved 

27 upon the macrogrid and control volumes using the initial and boundary 

28 conditions and properties dependent upon the volume fractions, i.e., average 

29 density and average viscosity, to arrive at local velocity conditions for the 

30 control volumes. The microgrid cells within the microgrid are then advected in 

31 response to the calculated local velocity conditions in the control volumes with 

32 voids and overlaps of the microcells occurring in the microgrid. The location 
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1 of each of the advected microcells is tracked using the unique enumeration 

2 previously assigned to each of the microgrid cells. Preferably, this unique 

3 enumeration includes the use of prime numbers, such as those generated by 

4 an Eulerian quadratic number generator. 
5 

6 Overlapped microgrid cells are then reallocated within the microgrid so that 

7 only one microgrid cell is located in each space of the microgrid to effectively 

8 conserve mass and satisfy local fluid fraction gradient values. The contents of 

9 the microgrid cells can then be used in additional time steps or iterations to 

1 0 determine new values for volume fractions of the fluid materials in the 

1 1 microgrid cells. The equations of fluid motion can then be iteratively solved 

12 again using updated properties dependent upon the volume fractions. This 

13 process of iteratively determining new fluid fractions in the control volumes, 

14 solving equations of fluid motion to arrive at new velocity fields, advecting 

15 microcells in accordance with the newly derived velocity fields, and then 

16 reallocating the microcells to comport with constraints such as conservation of 

1 7 momentum and conservation of mass is repeated over many time steps to 

1 8 simulate the fluid flow. 
19 

20 A method for tracking cells in a fluid dynamics computation is now described. 

21 Unique identifier numbers are assigned to microgrid cells located in a grid. 

22 The unique identifiers are associated with respective N fluid materials. N may 

23 be 2 or more materials. Preferably, the unique identifiers are prinie numbers, 

24 such as those generated by an Eulerian quadratic number generator. 
25 

26 The cells are advected within a grid in response to local velocity conditions 

27 such that some of the cells overlap one another while voids are also created 

28 in the grid. Grid locations with overlapped cells are identified by determining 

29 that a combination of the unique identifiers coexist in that grid location. The 

30 fluid materials of the microgrid cells disposed in a grid location can be 

31 determined using modular arithmetic when using unique identifiers which are 

32 prime numbers. The cells can then be reallocated to comport with specified 
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1 conservation of mass and momentum restrictions with one unique identifier 

2 again occupying each space in the microgrid. Volume fractions of fluid 

3 materials in specific areas of a computational space can then be computed to 

4 determine properties, such as average density and average viscosity, which 

5 are to be used in a next time step to again solve the time-dependent fluid 

6 equations of motion. 
7 

8 It is an object of the present invention to provide a method for simulating fluid 

9 flow that can easily handle two or more generally immiscible fluid materials 
1 0 and track the interfaces between the different materials. 

11 

12 It is yet another object to provide a method for tracking two or more materials 

13 and their associated interfaces based on the creation of elemental pieces or 

14 microgrid cells of fluid material and the use of a unique tagging system based 

15 on prime numbers, 
16 

17 It is still yet another object to provide a method for tracking N materials and 

18 their associated interfaces by uniquely identifying and tracking, with a single 

19 concentration function C, N fluid materials wherein N may be 2 or greater. 
20 

21 BRIEF DESCRIPTION OF THE DRAWINGS 

22 

23 These and other objects, features and advantages of the present invention 

24 will become better understood with regard to the following description, 

25 pending claims and accompanying drawings where: 
26 

27 FIG. 1 is a flowchart of the overall steps taken in modeling the flow of N fluid 

28 materials and tracking the interfaces between the fluid materials; 
29 

30 FIG. 2 illustrates a macrogrid including a macrocell or control volume which is 

31 overlayed with a 4 x 4 array of microgrid cells containing one of three distinct 

32 fluid materials which are identified by three corresponding prime numbers; 
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1 FIG. 3A shows a macrogrid including an array of 2 x 2 control volumes with an 

2 array of 3 x 3 microgrid cells overlaying each of the control volumes and FIG, 

3 3B is a corresponding chart including a matrix of prime numerals 

4 corresponding to the materials in the microgrid cells; 
5 

6 FIG. 4A shows an example of flux handling with prime enumerated microgrid 

7 cells showing a microgrid having overiapping cells and voids and FIG 4B is a 

8 corresponding chart including a matrix of prime numbers and products of 

9 prime numbers corresponding to the materials in the microgrid cells; 
10 

1 1 FIGS. 5A and 5B illustrate a visualization of a 2D scanning technique wherein 

12 FIG 5A shows a first scan and FIG. 58 shows a second scan which begins at 

13 a different starting point than the scan of FIG. 5A; 
14 

15 FIG. 6 is a flowchart describing steps for reallocating microgrid cells; 
16 

17 FIG. 7 is a flowchart describing steps in scanning a grid and reallocating 

18 microgrid cells into voids found in the microgrid; 
19 

20 FIG. 8 shows a definition for a theoretical flow in a rectangular cavity in which 

21 a plurality of fluid materials is allowed to swiri and mix; 
22 

23 FIG. 9 is a photographic representation of a test showing the swiri pattern in a 

24 rectangular cavity defined in FIG. 8; 
25 

26 FIGS. 10A and 10B show a shear driven recirculating cavity flow with three 

27 distinct numerical layers employing the microgrid algorithm of the present 

28 invention; 
29 

30 FIG. 11 shows a shear driven recirculating cavity flow with 3 distinct numerical 

31 layers employing the microgrid algorithm with a 2 x 2 refinement and a 

32 comparative 4x4 refinement at a time of 21 0 seconds; 
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1 FIG. 12 shows a shear driven cavity flow with 6 materials at time steps of 0, 

2 120, 240, and 360 seconds utilizing a 8 x 8 microgrid; and 
3 

4 FIG. 13 shows an artistic rendering of a shear driven cavity flow with 6 

5 colocated materials at 120 seconds. 
6 

7 BEST MODE FOR CARRYING OUT THE INVENTION 

8 

9 The present invention provides a method for simulating the flow of 

10 N-immiscibie fluid materials and tracking the position of the interfaces 

1 1 between the fluid materials. For the purpose of simplicity, an example utilizing 

12 three separate materials will be described. The method may easily be 

13 extended to model with a much larger number of fluid materials and 

14 interfaces. For example, using a particular number generator (Eulerian 

15 quadratic generator) and identification system for fluid materials, 39 materials 

16 may be used based on this prime number generator. Through use of other 

17 generators of unique identifiers, an even a higher number of materials and 

1 8 interfaces may be tracked. 
19 

20 An exemplary embodiment of a method for tracking N fluid materials and their 

21 interfaces, made in accordance with the present invention, is shown in a 

22 flowchart in FIG. 1. 
23 

24 A first step 1 10 is to create a macrogrid filled with macrocells or control 

25 volumes to define an overall computational space In which the flow of N fluid 

26 materials Is to be simulated. For example, fluid flow circulating In a 

27 rectangular cavity will be described in greater detail below. 
28 

29 Step 120 Includes applying initial and boundary conditions to the macrogrid 

30 and control volumes. In this exemplary embodiment, the boundary conditions 

31 include no-slip walls and a fixed tangential (or shear) velocity condition on the 

32 top surface of the cavity. Initial conditions are quiescent or zero velocity for 
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1 the fluids with a hydrostatic pressure distribution imposed in the vertical 

2 direction due to gravity forcing. Those skilled in the art of computational fluid 

3 dynamics modeling and simulation will appreciate that models can be made 

4 employing, by way of example and not limitation, other boundary and initial 

5 conditions such as symmetry planes, inflow and outflow boundaries, constant 

6 pressure surfaces, and an initial mean velocity and density distribution applied 

7 to the fluid materials. 
8 

9 A microgrid filled with microgrid cells is then overlaid upon the macrogrid and 

10 control volumes in step 130. FIG, 2 illustrates a single macrocell or control 

1 1 volume which is overlaid with a microgrid including a 4 x 4 array of microcells. 

12 In this particular preferred embodiment, all of the microgrid cells are of the 

13 same size and shape. The present invention could also be extended to use 

14 microgrid cells having different sizes and shapes. Microgrid cells offer a 

15 discrete and simple way to delineate materials in a computational space 

1 6 without the complexity of formal interface reconstoiction at every time step. In 

1 7 the microcell technique of the present exemplary embodiment of this 

18 invention, the resolution of an interface is predetermined by the size of the 

19 microgrid cells themselves and no cell by cell reconstruction is necessary. 
20 

21 N fluid materials are then assigned in step 140 to each of the microgrid cells. 

22 In this example, N= 3 as the flow of three fluid materials will be simulated. In 

23 FIG. 2, note that the microgrid cells are shaded to depict each of three 

24 different fluid materials 1 , 2 and 3. As a numerical representation of the 

25 physical system is created, a single prime number or identifier associated with 

26 a particular fluid material is assigned to each microgrid cell. For example, all 

27 microgrid cells that contain material N=1 are given the identifier 43. All 

28 microgrid cells that contain material N=2 are assigned the identifier 47 and 

29 cells containing material N=3 are labeled with the identifier 53. In this way, a 

30 prime representation of the system is built. Prime numbers less than 43 are 

31 not used in this exemplary embodiment because of a requirement that the 

32 square of any prime number be greater than the largest prime number 
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1 assigned to a material. This assignment renders the evaluated primes into a 

2 mathematically unique set which enhances the tracking of microcells. 
3 

4 Utilizing the same shading scheme as depicted in FIG. 2, FIG. 3 shows a 

5 macrogrid including a 2 x 2 array of control volumes each with a 3 x 3 array of 

6 microcells. Assume this is at an initial or first time step or t=0 seconds. The 

7 location of each material is accounted for in a microcellular fashion as shown 

8 in FIG. 3B utilizing a chart with a corresponding numerical representation of 

9 the microgrid cells and their assigned fluid materials. 
10 

1 1 A prime number generator is utilized in this preferred embodiment that will 

12 create enough unique primes to delineate every material in the simulation. 

13 This particular prime number generator has the following characteristics: 
14 

15 1 . the smallest prime number generated squared is greater than the 

16 largest prime number generated; and 

17 2. the generator creates a monotonically increasing set of primes. 
18 

19 By way of example and not limitation, a good 'all-purpose' generator that has 

20 these characteristics, and allows up to N=39 (where N is a positive integer), is 

21 the Eulerian quadratic generator, p(N) = + N + 41 . By the nature of how an 

22 algebraic generator works, some primes may be skipped. This particular 

23 generator creates a set of non-inclusive primes from 43 to 1601 as is 

24 demonstrated in Table 1 . Note that the largest prime generated, 1601 , is 

25 smaller than the smallest prime generated squared, 43^ = 1849. 
26 

27 Table 1 : Prime generator example. 
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In step 150, the volume fractions of the fluid materials in each of the control 
volumes are calculated. The prime marking technique is employed to resolve 
the fluid fraction of each material in each control volume. Fluid fraction, , 

may be recovered as a volume centered value at any time by using the 
formula: 



fN=QN 



V 

* MICRO 
V, ^MACRO J 



wherein is the number of microgrid cells occupied by material N, V^„cro 's 
the volume of a single microcell and V^acro 's the volume of a single control 
volume. 

Given the example of FIG. 2, and assuming evenly distributed microgrid cells 
in the 1 X 1 control volume. fi= 5/16 = 0.3125; h - 5/16= 0.3125 and fa = 6/16 
= 0.3750. 

Average volume centered values of density p and viscosity ^ , for example, 
in each control volume can then be recovered through: 

P = ^fuPu and // = 2^fN//N respectively. 

N N 

In step 160, equations of fluid motion are solved for each of the macrogrid 
control volumes to arrive at a velocity field for each of the control volumes (as 
well as other conserved variables such as pressure, turbulent kinetic energy, 
and temperature). The solution of these equations depends upon the volume 
centered values of density p and viscosity ju for each of the control volumes 
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1 (or, in general, all fluid property variables). These equations rely upon 

2 principles of conservation of mass, conservation of monnentum. and 

3 conservation of energy to derive velocity fields across each of the control 

4 volumes during each time step of the simulation. More particularly, the 

5 general equations of fluid motion include (written in tensor form): 

6 > 

7 5p/9t + 5pUj/5Xj = 0 (continuity equation) 
8 

9 p dujdl + pUj 5Ui/5Xj = - dP/dx\ + ay + pgi (momentum equation) 

10 

1 1 Cp dpT/dl + Cp 5pUiT/axi = dP/dt + Us aP/5Xi + 5/axi(k dT/dx{) + |iO + Q 

1 2 (energy equation) 
13 

14 where: Uj is the velocity is the i-coordinate direction, t is time, Xj is the spatial 

15 coordinate in the i-direction, P is pressure, ay is the strain rate tensor, Cp is 

16 constant pressure specific heat, T is temperature, k is thermal conductivity, <D 

17 is viscous dissipation, and Q is internal heat sources. 
18 

19 For simulation purposes, these equations are discretized into algebraic 

20 equations which are readily solved by iterative numerical techniques. 

21 Specifically, the discretized conservation equations for all control volumes are 

22 solved in an iterative manner. For example, when the x momentum equation 

23 is solved, all other unknowns are fixed (such as y and z momentum (in the 3D 

24 case), pressure, and density) to their most current iterative values, resulting in 

25 the only unknown variable being the x momentum. The resulting set of linear 

26 equations for the x momentum is then solved with a conjugate gradient type 

27 method for example. The x momentum field is then updated and this process 

28 is repeated for the y momentum, then z momentum, etc., working the iteration 

29 through each unknown variable. At the end of the iterative loop, all unknowns 

30 have effectively been updated to the next iteration level. This iterative cycle is 

31 continued until a set of convergence criteria are achieved (based on residuals 

32 and iterative value changes), and then the time step is complete and the 
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1 variable fields (velocity, pressure, temperature, etc.) are effectively updated to 

2 the next time and the above procedure is repeated again for the next time 

3 step. As a result of solving these equations, the velocity field (Vx, Vy, and Vz) in 

4 each control volume is determined as well as all other flow field variables 

5 (pressure, etc). 
6 

7 Those skilled in the art will appreciate that knowing the volume fractions, and 

8 subsequently average volume centered values of density p and viscosity // 

9 for each of the control volumes (as well as all other flow field variables), many 

10 different alternative solution schemes may be used to calculate the velocity 

1 1 fields or other values needed in order to calculate movement of the microgrid 

12 cells. The present invention provides the advantage that numerous different 

13 fluids materials may be used in a simulation with the volume fractions of the 

14 control volumes being readily calculated. From these volume fractions, 

15 average parameters needed to solve equations of fluid motion can be readily 

16 derived. Consequently, a wide number of methods or algorithms can be used 

17 to estimate velocity fields or other variables which are useful in computational 

18 fluid dynamic simulations. These parameters would then be used to advect or 

1 9 move microgrid cells as will be described below. 
20 

21 Step 170 includes advecting the microgrid cells within the macrogrid and each 

22 of the control volumes in accordance with the iterative velocity fields 

23 calculated in step 160. This displacement of microgrid cells will result in some 

24 of the microgrid cells overlapping one another in some of the microgrid 

25 locations and other microgrid locations being left void. FIG. 4 illustrates such 

26 a microgrid containing overlapping and void microgrid locations. 
27 

28 The microgrid cells move only as a complete unit and individually represent a 

29 single material. A fundamental premise of the microgrid scheme is that a 

30 material may not become less than 1 microgrid in size. The microgrid is the 

31 fundamental length or volume scale in the simulation. Each microgrid cell 

32 reacts to its associated local control volume or macrogrid velocity conditions 
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1 and may subsequently be displaced to another control volume after 

2 movement is rectified. 
3 

4 Movement or displacement of a microgrid cell in a particular direction is 

5 determined by the velocity field of the control volume in a particular direction 

6 (x, y or z) and the duration of movement, i.e., the length of the time step. For 

7 example, if the average velocity of a control volume in the x-direction, as 

8 determined in step 160, is 0.5 meters/second and the duration of a time step 

9 is 1 second, then each of the microgrid cells in that control volume would want 

10 to be moved 0.5 meters along the x-direction. If the size of each microgrid 

1 1 cell is 0.5 meters, then each microgrid cell in that particular control volume 

12 would be displaced by one microgrid cell. Similarly, the microgrid cells in the 

1 3 control volumes may also move in the y- and z- directions, if a 3 dimensional 

14 model were used. 
15 

1 6 Movement of the microgrid cells is only allowed if the microgrid cell has 

1 7 reached the threshold of moving at least one microgrid cell unit. In this 

18 exemplary embodiment, when a microgrid cell experiences a velocity 

19 condition too small to warrant movement, this information is flagged and the 

20 next iteration step is performed and so on until movement is imminent based 

21 on the local velocity field. If the velocity field is not of sufficient strength to 

22 warrant movement of a microgrid cell, then it simply remains in its current 

23 location. In this way, an accurate accounting of the effects of velocity can be 

24 maintained while utilizing a somewhat coarse microgrid. 
25 

26 As the advection step begins, an interim array is used to receive advected 

27 microgrid cell information. This interim array is initialized to values of "1", and 

28 the size of the interim array corresponds to the size of the microgrid cellular 

29 field. During the advection step, new microgrid cell values are constructed via 

30 a multiplicative association of primes and then stored in the interim array at 

31 the appropriate locations. For example, a microgrid cell receiving one copy of 

32 a single material (such as material 1 ) would be represented by 1 x 43 = 43. A 

-14- 



1 microgrid cell receiving 2 different materials may be numerically represented 

2 as 1 X 43 X 53 = 2279, which indicates the presence of both materials 1 and 3 

3 in the microgrid cell. Similarly, a cell may receive two copies of the same 

4 material (e.g. 1 x 43 x 43 = 1849). If no material is advected into an interim 

5 array location, then the value remains as "1" and indicates that void is 

6 present. Again, FIG. 4A illustrates the spatial distribution of prime numbers in 

7 a macrogrid after the advection step 170 is complete. (Note that the original 

8 distribution of microcell fluid materials before the advection step is shown in 

9 FIG. 3.) Where there is a void, the array contains a numeral "1". Where there 

10 is an overlap of microgrid cells, the product of the material identifers (primes) 

11 of the overiapping microgrid cells is given. 
12 

1 3 After every iteration within a time step the microgrid cells must be reallocated, 

14 step 180, so that there are no remaining overiapping cells or voids in the 

1 5 microgrid. Microgrid locations that contain more than one identifying prime 

1 6 number or material must be stripped to a single prime number with the other 

17 identifiers shifted to new locations. Microgrid locations void of any prime 

18 identifier must be filled. The reallocation scheme preferably is designed to 

1 9 achieve the objectives of insuring conservation of variable fields (mass, 

20 momentum, and energy). 
21 

22 A correction algorithm is designed to detect values that are larger than the 

23 largest prime in the system, which are representative of overiaps of microgrid 

24 cells. These few numbers are then evaluated to ascertain the contents of the 

25 associated microgrid cells. By controlling what primes are available to create 

26 interim microgrid cell values, the evaluation of overiapping and void microgrid 

27 spaces is reduced to a single greater than or less than comparison. 

28 Consequently, when a large number is detected, the possible factors of that 

29 number are contained in a small subset of the known primes. Specifically, 

30 these are the primes created by the Eulerian quadratic generator, p(N), 

31 defined in Table 1 . All together, this technique greatly reduces the 
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1 computational cost to evaluate the contents and reorganization of microgrid 

2 cells. 
3 

4 For example, given the numeric representation in FIG. 4B, the algorithm 

5 detects that the number 2279 is larger than the largest prime in the system 

6 (which in this case is 53) and begins a scheme that delineates the contents of 

7 the cell, breaking the integer 2279 into its prime components that represent 

8 the materials therein. Specifically, the algorithm exploits the concept of 

9 modular arithmetic until a solution of 0 is calculated. Recall from the rules of 
10 modular arithmetic that (2279 mod 53) = 43. And, (43 mod 43) = 0. 

11 

12 A preferred exemplary reallocation procedure is described in Flowchart 2 of 

13 FIG. 6. Again, those skilled in the art will realize that many other alternative 

14 reallocation procedures or algorithm could be used to rearrange the 

15 overlapping microgrid cells and voids to establish a microgrid having a single 

16 fluid material or microcell in each microgrid location, as was the .case in the 

17 initial microgrid cell arrangement shown in FIG. 3. The different reallocation 

18 procedures can be used to achieve different reallocation objectives such as 

19 using subiterations to balance material distributions. 
20 

21 The first step 220 in FIG. 6 is to scan the entire microgrid field for "1" *s or 

22 void locations. Recall that the temporary grid that holds the transported or 

23 fluxed identifiers (primes) is initialized with 1's. If no primes flux into a 

24 microgrid location, the value "1" is maintained. These locations correspond to 

25 a microcellular void and a sum of their number gives a maximum number of 

26 overlap sites that subsequently sizes a caching array. It is possible that a 

27 single location may hold more than one overlap so void locations may 

28 outnumber overlapped locations. The caching array is used to record excess 

29 or stripped identifiers from microgrid cells during the reallocation process. 
30 

31 Next in step 230, the microgrid is scanned for microgrid cells that contain 

32 overlapped single materials. These are locations that contain multiple copies 
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1 of the same prime number. In this situation, the excess primes are stripped 

2 (step 240) from the microgrid cells using modular arithmetic and then placed 

3 into the caching array, while one copy is left behind to fill the location. 

4 Mathematically, the following evaluation occurs: 
5 

6 If: / > p{N^ ) and MOD{f, p{n)) = 0 where 

7 / > p{N), p(N) = A^' + AT + 41, and /, E Z\ 

8 N = {N,,N,,N,,,..,N^J. 



9 then there are multiple copies of the same prime (or material). 

10 In step 250, subsequently the "1" 's locations are quickly evaluated against 

1 1 their gradient condition. The gradient condition is evaluated by solving VfN = 

12 0 . The bulk void is filled using available primes from the caching array by 

13 evaluating the gradient condition (in 2D or 3D space). This means that If all 

14 the microgrid cells in closest proximity (i.e. every touching microgrid cell) are 

15 the exact same prime number, and there is an equal prime number available 

16 in the caching array, the "1" is immediately replaced with the corresponding 

17 prime. This reduces the number of locations that must be evaluated in a 

18 conservation of mass routine, to be described below. 
19 

20 At this point the microgrid system may still contain overlaps that hold different 

21 primes as well as still having void or "1" locations. The "V locations typically 

22 develop at interfaces, grid edges, or in areas of large scale (relative to the 

23 mesh) bulk motion. The overlaps are all interfacial areas. 
24 

25 Once again, the microgrid is scanned for cells that contain more than one 

26 prime (or different materials). This is done in step 260 by exploiting the 

27 property of the Eulerian generator, p(N) = + N + 41 (N = 1,2, ..39), such that 

28 the smallest prime squared is larger than the largest prime. 

29 If: / > p{N^ ) and MOD{/,p{n))>0 where 

30 / > p(N), p(N) = iV' + + 41, and /,A^ G Z^ 

31 N = {N,,N,,N,,,..,N^l then. 



-17- 



1 

2 these locations correspond to microgrid cells with interfacial overlap. 
3 

4 When these cells are located, all of the prime identifiers are removed from the 

5 location. An evaluation is made as to which prime most appropriately fills the 

6 location. This is done by an evaluation of local fluid fraction gradient and a 

7 momentum value assigned to each material based on the local velocity 

8 condition. More particularly, the overlapped cells are evaluated based on 

9 gradient, V/^ and momentum, pf^\v\, criteria, such that a physical motivation 

10 is established to identify the correct material distribution. The most ideal 

1 1 candidates for primes are chosen and the excess primes are stripped to the 

12 caching array. 
13 

14 Once all of the overlaps are corrected, all that is left are microgrid cellular 

15 locations that contain a "1" and the remaining material identifiers in the 

16 caching array. In step, 280, the remaining voids are evaluated using 

17 maximum material gradient condition, momentum considerations and material 

18 content of the caching array. Ideally, a particular material or prime is selected 

19 such that its reallocation provides a maximum or optimal gradient condition. 

20 Mathematically, this means that Pn, such that min (Vfw) for N = 1 , 2, 3 This 

21 choice of material is preferably overridden if a microcell exceeds a 

22 predetermined momentum criteria, i.e. momentum = mass of microcell x 

23 velocity which is greater than a predetermined threshold for momentum. If so, 

24 the high momentum microgrid cell(s) are located in the void. In this way, 

25 individual high momentum microcells are allowed to "pierce" otherwise 

26 continuous flows of other materials. 
27 

28 FIG. 7 shows a third flowchart which summarizes step 280 used in evaluating 

29 and filling these remaining voids. In step 320, the microgrid voids are 

30 evaluated for maximum gradient condition. More specifically, the maximum 

31 gradient condition is evaluated by scanning linearly for locally dissimilar 

32 material as described in step 330a. FIG. 5A illustrates that from a particular 

-18- 



1 starting "S" void location, neighboring microgrid cells are scanned in 8 linear 

2 directions. The scan begins in the direction with the greatest gradient. If the 

3 prime selected is available in the caching array, a replacement is made in the 

4 particular void. 
5 

6 In step 340a, the needed prime is extracted from its closest neighbor and is 

7 replaced with a prime from the caching array. In step, 350a the void is 

8 replaced with the just extracted prime. 
9 

10 If a desired prime is not available in the caching array, a scan, step 330b, is 

1 1 made starting from the particular microgrid void or starting position "S" and 

12 moving out in every linear direction (8 in 2D, 26 in 3D) where an attempt is 

13 made to find the prime locally. Once found, the microgrid is adjusted in that 

14 direction to compensate for a material imbalanced distribution. In step 340b, 

15 the needed prime from the closest neighbor is extracted and a void is created 

16 in a new location. The original void is then replaced with the extracted prime, 

17 step 350b. The 'new' void is then evaluated for maximum gradient condition 

18 in step 360b. In step 370b, the grid is scanned linearly from a new void area 

19 for local dissimilar material that already exists in the caching array. Steps 

20 330a, 330b and 330c are then repeated. 
21 

22 Step 330c provides a procedure to address the issue of an unsuccessful grid 

23 scan for dissimilar materials. Using a randomizer, the starting location for the 

24 scan is moved in step 340c. In step 350c, the contents of microgrid cells are 

25 exchanged. Again the steps may be taken as described in steps 330a, 330b 

26 or 330c. 
27 

28 Usually what results is a slight local reorganization that conserves mass. If 

29 this is not possible, a mass error is accrued and through the progress of a 

30 calculation, the total mass error is determined and used as a gauge to 

31 evaluate the overall simulation accuracy. 
32 
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1 Once broken down, a separate algorithm, described in Flowchart 3 of FIG. 6 

2 evaluates the local condition Vf^ for each N and designates the best 

3 candidate prime for the previously overiapped microgrid cell. The prime NOT 

4 replaced in the once mixed cell is then redistributed to a next best location. 

5 Ideally, this location is the closest cell in the direction of its strongest material 

6 gradient. The algorithm is designed to scan in this direction in search of an 

7 empty cell that will be best filled by the prime. If no such cell can be found, in 

8 2D, 8 directions are scanned as shown in FIG. 5A. If still no candidate can be 

9 located, the initial search is abandoned. 
10 

11 A small adjustment to the starting location of the search is made and the 

12 multi-directional scan is resumed (see FIG. 5B). The user can define how 

13 many times these small adjustments will be tolerated before an "interfacial 

14 bump", for lack of a better description, is implemented. What happens is this: 

1 5 the direction associated with the largest material gradient is rescanned for a 

16 closest interface. At the interface, a material exchange is made, replacing the 

17 closest interi'acial cell with the local material at hand and then using this 

18 location as a starting point to seek out a new site for the material that was 

19 displaced, using the original searching algorithm. After several iterations of 

20 this sequence, all empty microgrid cells and overiaps are corrected. 
21 

22 If sufficient iteration steps have occurred to achieve convergence, the 

23 iteration is ended. If additional time steps are required, then step 150 is 

24 repeated using the reallocated cells and their primes from step 1 80 to 

25 redetermine the volume fractions in each of the control volumes. Steps 160- 

26 190 are repeated until it is determined that sufficient time steps have occurred 

27 and the simulation is ended in step 200. 
28 

29 Computational Example 1 . 
30 

31 Circulating rectangular cavity flow offers a rigorously investigated, 

32 complicated flow pattern for the reader's review. In this example, the cavity is 

-20- 



1 initialized with a stratified pattern of three similar nnaterlals, differing only in 

2 prime designation not material properties. In this manner the similarity of the 

3 flow form is qualified by literature comparison to a single fluid flow, as well as 

4 to illustrate the use of three numerically distinct materials in the microgrid 

5 algorithm. 
6 

7 Rhee et. al. (1984) and Freitas (1986) describe a similar recirculating cavity 

8 flow such that given the schematic of FIG. 8, a plate drives a continuous 

9 shear motion from front to back at the top of the cavity. FIG. 8 also illustrates 

1 0 the formation of secondary eddies in three of the four corners of the cavity 

1 1 along a center plane. FIG. 9 is a photograph of Rhee's experimental flow 

12 along a similar plane for comparison. FIG. 10 is a rendering of the microgrid 

13 algorithm output for the three stratified materials described in the previous 

14 paragraph. 
15 

16 As shown in FIG. 10, clearly a portion of the bottom layer is trapped in the 

1 7 lower two eddies, while the balance of the material travels through the 

18 recirculation pattern. The middle layer is most heavily mixed as it is quickly 

19 taken up in the primary circulation cell then partially fragmented as it makes 

20 its first full rotation. The top layer is partially trapped in the upper left hand 

21 corner eddy with the balance following the primary circulation cell. It is 

22 obvious that the microgrid algorithm is correctly capturing the proper features 

23 of the recirculating cavity flow physically documented in FIG. 9 by Rhee et. al 

24 (1984) and computationally rendered by Freitas (1986), while simultaneously 

25 tracking three numerically separate materials. 
26 

27 To capture this refined detail, the simulation employs a fairly coarse mesh. 

28 The physical dimension of the cavity portrayed in FIG. 10 is 150 mm x 150 

29 mm. The computational macrogrid is 126 x 126. The microgrid cell mesh is 4 

30 X 4. Numerical experiments verify that an increase in microgrid mesh while 

31 maintaining the original computational mesh yields more refined results with a 

32 minimal increase in computational expense. 
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1 

2 EXAMPLE 2 
3 

4 As an example, this same simulation is offered with 126 x 126 computational 

5 mesh and a microgrid cell mesh of 2 x 2 (1/4 of the microgrid resolution 

6 shown in the above simulation of FIG. 10) in a side by side comparison at 

7 time step =210. This time step is chosen because it shows adequate detail of 

8 all the eddy features. 
9 

10 Because the algorithm employs an additive and collective response measure 

11 to the effect of velocity over multiple time steps, the more coarse the microgrid 

12 cell refinement is, the less immediate physical response one can see in areas 

13 of relatively low velocity. Similarly, because of the block-like movement of 

14 microgrid cells, to capture areas of high thinning, microgrid cell size must be 

15 reduced. This Is most evident in the middle layer, where material should form 

16 thin trails in response to the shear motion of the primary circulation cell. Also, 

1 7 while the overall response of the system at coarse micromesh is good, it may 

18 be noted that interfacial smoothness and subtle feature resolution Is increased 

19 with the refinement of the microgrid. It is important, however, to realize that at 

20 no time is the computational mesh adjusted. These calculations are coupled 

21 to the microgrid via a fluid fraction value only, and the computational mesh is 

22 completely independent of microgrid cell size. The best computational value, 

23 therefore, utilizes a coarse computational mesh with a refined microgrid mesh. 
24 

25 As the only algorithmic limitation of the microgrid process demonstrated here 

26 with p(N) = + N + 41 is a numerical cap of 39 materials, it can easily be 

27 demonstrate that the simulation methodology can be used with a larger 

28 number of materials. Of course, the computational grid selection must be 

29 chosen with prudence. Just as in a 2 material simulation, it is not wise to 

30 represent a single material over too few computational cells. 
31 

32 
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1 EXAMPLE 3 
2 

3 The results shown in FIG. 12 utilize a 126 x 126 computational mesh with an 

4 8x8 microgrid mesh. FIG. 13 is an artistic rendering of the 6 materials in the 

5 same view at time step =120, 
6 

7 While in the foregoing specification this invention has been described in 

8 relation to certain preferred embodiments thereof, and many details have 

9 been set forth for purposes of illustration, it will be apparent to those skilled in 

10 the art that the invention is susceptible to alteration and that certain other 

1 1 details described herein can vary considerably without departing from the 

12 basic principles of the invention. 
13 
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